function IRF = BVARir(VAR, sigma_draw, beta_draw, nsteps)
% =======================================================================
% Compute the Choleski IRFs for a VAR model with a given VCV matrix and a
% given matrix of estimated parameters. ONLY THE ONE shock specified in
% impulse is computed. This program can be used to get the IR of a VAR 
% estimated with Bayesian methods (with BVAR_drawpost.m for example). 
% =======================================================================
% IRF = BVARir(VAR, sigma_draw, beta_draw, nsteps)
% -----------------------------------------------------------------------
% INPUT
%   VAR        : structure, result of VARmodel function
%   sigma_draw : draw from the posterior distribution for sigma
%   beta_draw  : draw from the posterior distribution for beta
%   nsteps    : number of nsteps to compute the IRFs
%
% OUPUT
%   IRF(:,:,:) : matrix with nsteps, variable, shock
% =======================================================================
% Ambrogio Cesa Bianchi, May 2012
% ambrogio.cesabianchi@gmail.com



%% Retrieve parameters and preallocate variables
%===============================================
c_case   = VAR.c_case;
nvars    = VAR.neqs;
nlags    = VAR.nlag;
low_chol = chol(sigma_draw)';
IRF      = zeros(nsteps,nvars,nvars);

%% Compute the impulse response
%===============================

for mm=1:nvars

    % check for constant, trend, and/or exogenous variables and remove them
    switch c_case
        case 0
            beta_draw_t = beta_draw';
            big_gamma = [beta_draw_t(:,1:nvars*nlags)
                eye(nvars*(nlags-1)) zeros(nvars*(nlags-1),nvars)];

        case 1
            beta_draw_t = beta_draw';
            big_gamma = [beta_draw_t(:,2:nvars*nlags+1)
                eye(nvars*(nlags-1)) zeros(nvars*(nlags-1),nvars)];

        case 2
            beta_draw_t = beta_draw';
            big_gamma = [beta_draw_t(:,3:nvars*nlags+2)
                eye(nvars*(nlags-1)) zeros(nvars*(nlags-1),nvars)];
    end

    big_gamma_pot = eye(size(big_gamma,1));
    
    % Initialize the impulse response vector
    imp_res = zeros(nvars, nsteps);
    
    % Create the impulse vector
    impulse = zeros(nvars,1); 
    impulse(mm,1) = 1;                 % 1 stdev shock
%     impulse(mm,1) = 1/low_chol(mm,mm); % 1 unity shock

    % First period impulse response (=impulse vector)
    imp_res(:,1) = low_chol*impulse;
    big_impulse  = [(low_chol*impulse)' zeros(1, nvars*(nlags-1))]';
    
    % Recursive computation of impulse response
    for kk = 2:nsteps
        big_gamma_pot = big_gamma_pot * big_gamma;
        big_imp_res   = big_gamma_pot * big_impulse;
        imp_res(:,kk) = big_imp_res(1:nvars);
    end
    IRF(:,:,mm) = imp_res';
end
